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Abstract. A pathway-based mean-field theory (PBMFT) was recently pro- 
posed for E. coli chemotaxis in [G. Si, T. Wu, Q. Quyang and Y. Tu, Phys. 
Rev. Lett., 109 (2012), 048101]. In this paper, we derived a new moment sys- 
tem of PBMFT by using the moment closure technique in kinetic theory under 
the assumption that the methylation level is locally concentrated. The new 
system is hyperbolic with linear convection terms. Under certain assumptions, 
the new system can recover the original model. Especially the assumption on 
the methylation difference made there can be understood explicitly in this new 
moment system. We obtain the Keller-Segel limit by taking into account the 
different physical time scales of tumbling, adaptation and the experimental 
observations. We also present numerical evidence to show the quantitative 
agreement of the moment system with the individual based E. coli chemotaxis 
simulator. 



1. INTRODUCTION 

The locomotion of Escherichia coli {E. coli) presents a tumble-and-run pattern 
([4]), which can be viewed as a biased random walk process. In the presence of 
chemoeffector with nonzero gradients, the suppression of direction change (tumble) 
leads to chemotaxis toward the high concentration of chemoattractants ([ijjs]). 
Great efforts have been put into understanding the chemotactic sensory system 



of E. coli ([15 30 32 ). The chemotaxis signaling pathway belongs to the class of 
two-component sensory system, which consists of sensors and response regulators. 
The chemotaxis sensor complex is composed of transmembrane chemo-receptors, 
the adaptor protein CheW, and the histidine kinase CheA. The response regulator 
CheY controls the tumbling frequency of the flagellar motor (|16 ). Adaptation 
is carried out by the two enzymes, CheR and CheB, which control the kinase 
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activity by modulating the methylation level of receptors ([32 ). Because of slow 
adaptation process, receptor methylation level serves as the memory of cells, and 
cells decide whether to run or tumble by comparing receptor methylation level to 
local environments. 

In the modeling literature, bacterial chemotaxis has been described by the Keller- 
Segel (K-S) model at the population level ([20]), where the drift velocity is given 
by empirical functions of chemoeffector gradient. It has successfully explained 
the chemotaxis phenomenon in slowly changing environments ( [29| ), however fails 
to make good predictions in rapidly changing ones ([34]) and the volcano effects 
([9||26]). Besides that, the K-S model has also been mathematically proved to 
present nonphysical blowups in high dimension when initial total mass reaches the 
critical level (|6]-[8]). In order to understand bacterial behavior from the individual 
dynamics, kinetic models have been also developed by considering the velocity-jump 
process ( [3{ [T7|[28) ) , and the K-S model can be systematically derived by taking the 
hydrodynamic limit of kinetic models {e.g. 10 M]). All the above mentioned mod- 
els are phenomenological and do not take into account the signaling transduction 
and adaptation process. 

Nowadays, modern experimental technologies have been able to quantitatively 
measure the dynamics of signaling pathways of E. coli ( [2j[ll|[24j[27j ) , which has led 
to successful modeling of the pathway dynamics ( [2l|[23|[3l] ) . These works make 
possible the verification of predictive agent-based models that include the intracel- 
lular signaling pathway dynamics. It is of great biological interest to understand 
the molecular origins of chemotaxis behavior of E. coli by deriving population-level 
model based on the underlying signaling pathway dynamics (|12[|25|). Particu- 



larly in 25 , the authors developed a pathway-based mean field theory (PBMFT) 
that incorporated the most recent quantitatively measured signaling pathway, and 
explained a counter-intuitional experimental observation which showed that in a 
spatial-temporal fast-varying environment, there exists a phase shift between the 



dynamics of ligand concentration and center of mass of the cells 34 . Especially, 



when the oscillating frequency is comparable to the adaptation rate of E. coli, the 
phase shift becomes significant. Apparently this is a phenomenon that can not be 
explained by the K-S model. 

In this paper, we study PBMFT for E. coli chemotaxis from a mathematical 
point of view. Specifically we derive a new moment system of PBMFT using the 
moment closure technique in kinetic theory. The new system is hyperbolic with 
linear convection terms. Under certain assumptions, the derived moment system 
gets to the original model in 25 , and especially the assumption on the methy- 
lation difference made in 25 can be understood explicitly in this new system. 
Taking into account the different physical time scales of the tumbling, adaptation 
and experimental observation, we connect the moment system to the K-S model 
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(in the parabohc scaUng). The agreement of the moment system with the signal- 
ing pathway-based E. coli chemotaxis agent-based simulator (SPECS [l9]) will be 
provided numerically in the environment of spatial-temporal varying ligand con- 
centration. 

The rest of the paper is organized as follows. We introduce the pathway-based 
kinetic model incorporating the intracellular adaptation dynamics in Section [2j In 
Section [3] assuming the methylation level is locally concentrated, we are able to 
build the moment system by using the moment closure technique in one dimen- 
sion. Furthermore, the modeling assumption will be justified both analytically and 
numerically. Section |4] illustrates why K-S model is valid in the slow varying en- 
vironments. We also give the connection of the moment system to the PBMFT 
model proposed in 25 , and provide the quantitative agreement of the moment 
system with SPECS numerically. Two-dimensional moment system is derived in 
Section [5j and we make conclusive remarks in Section |6] 

2. Description of the kinetic model 
We shall start from the same kinetic model used in 25j, which incorporates the 



most recent progresses on the chemo-sensory system ([24,31 ). The model is a 
one-dimensional two-flux model given by 

dP+ _ d{voP+) d{f{a)P+) z{m) 
^ > dt ~ dx dm 2 ^ 

(2.2) dP- ^ d{v^P~) d{f{a)P-) ^ z{m) 

dt dx dm 2 

In this model, each single cell of E. coli moves either in the or "— " direction 
with a constant velocity vq. P^{t,x,m) is the probability density function for the 
cells moving in the "±" direction, at time t, position x and methylation level m. 
The intracellular adaptation dynamics is described by 

(2.3) ^ = /(«) = fcfl(l-a/«o), 

at 

where the receptor activity a(m, [L]) depends on the intracellular methylation level 
m as well as the extracellular chemoattractant concentration [L] , which is given by 

(2.4) {l + exp{NE)y\ 
According to the two-state model in ^l][23] , the free energy is 

1 + [L]/Ki \ 
1 + [L]/Ka)' 

In ( |2.3[ ), kji is the methylation rate, ao is the receptor preferred activity that satisfies 
/(ao) = 0, f'{ao) < 0. TV, mo, Kj, Ka represent the number of tightly coupled 
receptors, basic methylation level, and dissociation constant for inactive receptors 
and active receptors respectively. 



(2.5) E = -a{m-mo) + h{[L\), with /o([i]) = In 



4 



GUANGWEI SI, MIN TANG, AND XU YANG 



We take the tumbling rate function z(m, [L]) in (25] , 

(2.6) z = zo + T-\a/ao)", 

where zq, H, t represent the rotational diffusion, the Hill coefficient of flagellar 
motor's response curve and the average run time respectively. We refer the read- 
ers to [25| and the references therein for the detailed physical meanings of these 
parameters. 

More generally, the kinetic model incorporating chemo-sensory system is given 
as below, 

(2.7) dtP = -V ■ V^P - drn{fia)P) + Q{P, z), 

where P(i, a;, v, m) is the probability density function of bacteria at time i, position 
X, moving at velocity v and methylation level m. 

The tumbling term z) is 
(2.8) 

Q{P,z)— -j- z{m,[L],v,v')P{t, x,v' ,m) dv' — + z{m,[L],v' ,v) dv'P{t,x,v,m), 
Jn Jn 

where n represents the velocity space and the integral 



1 



where = / dv, 
o Jn 



denotes the average over fl. z{m, [L], v, v') is the tumbling frequency from v' to v, 



which is also related to the activity a as in (2.6). The first term on right-hand side 



of (2.8 1 is a gain term, and the second is a loss term. 



3. One-dimensional mean-field model 



In this section, we derive a new moment system of PBMFT from (2.1)-(2.2) 
based on the the assumption that the methylation level is locally concentrated. 
This assumption will be justified by the numerical simulations using SPECS and 
the formal analysis in the limit of /c^ — > oo. To simplify notations, we denote /q^°° 
by / in the rest of this paper. 

3.1. Derivation of a new moment system of PBMFT. Firstly, we define the 
macroscopic quantities, density, density flux, momentum (on m) and momentum 
flux as follows. 



(3.1) p{x,t)= j {P+ +p-)dm, Jp{x,t) = j vo{P+ - P-)dm; 

(3.2) q{x,t) = J m{P+ + p-)dm, Jq{x,t) J vom{P+ ~ P') dm. 
The average methylation level M{t,x) is defined as 
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(3.3) 



M 



For simplicity, we also introduce the following notations 

dZ dz 

Z = z[M{t,x)), 

(3.4) 



Z = z{M{t,x)), 
F = f(a{Mit,x))y 



dm dm 
OF _ df 
dm dm 



M 



m=M 



Assumption A. We need the following condition to close the moment system, 

j{m-MYP^ dm 

/P± dm ^ ^' 

Remark. Physically this assumption means, the distribution functions is lo- 
calized in m, and the variation of averaged methylation is small in both moving 
directions "±". 



Integrating (2.1)+(2.2| with respect to m yields the equation for density 

= 0. 



dp dJp 



dt dx 



Integrating uox(2.1|— vqX.(2.2) with respect to m produces 



5Jp 
dt 



-v^— Uo / z(m)(P^ — P ) dm 

dx 



2 dp 
"'dx 



vo I [ z{M) + 
dZ 



dm 



{m~M) \{P+ - p-)<lm 



where we have used Assumption [A| in the second step and the notations in (3.3), 
(3.4) in the third step. 



Similarly, integrating mx(2.1)+ mx(2.2) with respect to m gives 

dq 
di 



dq dJq 



dx 

dJg 
dx 



+ / f{a){P+ + p-)dm 



+ / /(a)|m=M + 



dl 
dm 



=M 



{m- M) ]{P+ +p-)dr 



dx dm 

dJa 



dx 



Fp, 



where we have used an integration by parts in the first step and the definition of 
M in (3.3) in the last step. 
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Integrating wottix (2.1 )— vomx(2.2) with respect to m yields 

{m- M)]{P+ - p-)dm 



^ = - «o|^ + «o / /(a)(P+ - P-) dm-vo I z{m)m{P 



Wo / iz{m)m)\rn=M + 



d{z{m)m) 



dm 



dm 
OF 
dm 



{m- M)j{P+ - p-)dm 

M + (J, - Af Jp) 



m=M 

^ dx ' "'' ' dm^""^ "'^ dm 

where we have used Assumption [A] in the last step. 
Altogether, we obtain a closed moment system for p, Jp, q and Jq 



(Jg - MJp) - ZJq - ^M(J, - MJp), 



(3.5) 
(3.6) 
(3.7) 
(3.8) 



dp 


_dJp 


dt 


dx ' 


dJp 




dt 




dq 




dt " 


dx 


dJq 




dt 





— ZJp — — — 
dm 



(Jq^MJp), 



Fp, 



FJo 



dp d7 

— (J, - MJp) ZJq g^M{Jq - MJp). 



Remark. The Taylor expansion in m gives a systematical way of constructing high 
order moment systems. Please see the Appendix for the derivation of second-order 
moment system. 

3.2. Numerical Justification of Assumption |A] by SPECS. To justify the 
Assumption |A] we simulate the distribution of m with SPECS in an exponen- 
tial gradient ligand environment [L] ~ [i]o exp(Ga;). SPECS is a well developed 
agent-based E. coli simulator that incorporates the physically measured signaling 



pathways and parameters. We refer the readers to 19 for its detailed description. 
In the simulation, cells exiting at one side of the boundary will enter from the other 
side, and the methylation level is reset randomly following the local distribution 
of m at the boundaries. Under this boundary condition, the system will reach the 
steady state after a period of transient process. The steady state distributions are 
shown in Figure [l] In each of the subfigures, the horizontal and vertical axes repre- 
sent the position and the methylation level respectively. As shown in Figure [T] the 
distribution of methylation level is localized, and becomes wider when G increases. 
M± = / mP± dm are the average methylation levels for the right and left moving 
cells. One can also observe that M+ < M~ in the exponential increasing ligand 
concentration environment. This can be understood intuitively by noticing that 
the up gradient cells with lower methylation level come from left while the down 
gradient cells with higher methylation level come from right. 
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Figure 1 . The distribution of cells' receptor methylation level for 
G = 0.0005/im-i (a) and G ^ Qml^^lm'^ (b). The red dots 
represent cells moving to right while the blue ones represent those 
moving to left. are the average methylation levels for the 

right and left moving cells respectively. In the simulation, we take 
[i]o = bKj. Other parameters used in the SPECS are the same as 



those proposed in 19 . 



In the exponential environment, the numerical variations of m are almost uniform 
in X. The maximum of the methylation level variation in the simulation domain is 
defined by 



' /(m - M{x)f{P+ + P-)dm 



V J{P+ + P-)dm 

Assumption [X] is equivalent to the condition ct <C 1. As shown in Figure [2] a 
increases in G and decreases in kfj, but it is always small in the parameter regime 
we are interested in, i.e. Assumption [K\ holds in these cases. 

3.3. The localization of P^ in m in the limit of kji ^ 1. We show by formal 
analysis that the assumption J{m — M^P^ dm ^ 1 is true when the adaptation 
rate A:/} ;» 1. Denote 

(3.9) kR = 1/r,, /(a) - f^{a)M, 



then (2.1|-(2.2| become 



niO^ aKP+) i a(/,(a)P+) z + 

^^■^"^ dt " dx dm 2^^ ^ >' 

(3.11) 0F_^9J^_18if^^^ 
at ax rj am 2 
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Figure 2. The maximum variances a oi m for different G and 
kji. a increases in G for a given fcjj (a) and decreases in /c/j with 
fixed G (b), but they are aU smah in the parameter regime we are 
interested in. 



Integrating the above two equations with respect to m produces, for P^(t,x) = 
/cT ^^(^' ™) "^"^ (-^ arbitrary positive constant), 

r-R 

P^) dm 



(3.12) 



(3.13) 



dt 



9PR 
dt 



dx 



1 



2 



^-^(a(i?))P+(t,a:,i?) + -/^ (a(0))P+(t, a:, 0), 



dx 



P ) dm 



- -/^ (a(i?)) P- (i, X, P) + -/^ (a(0)) P- (t, J, 0) . 

The probability density functions satisfy P^{t,x,'m) > 0, Vm > 0, and thus 
P^(t,x) increases with R. 
We consider the regime 

(3.14) r?<l, and /^(a) - 0(1). 
Then when < 1, ( |3.12| )-( |3.13[ ) indicate for P e (0, +oo), 

(3.15) /„(a(P))P±(<,a;,P) -/„(a(0))P±(t,a:,0)+O(r,). 

We show by contradiction that when 77 — > 0, the boundary condition at m = has 
to satisfy /^(a(0))P='=(t, x, 0) 0, V(a;,t) e M x (0, +00). Otherwise, assume that 

(3.16) fn{a{0))P^{t,x,0) ^ G{t,x) ^0, for some (x, i) e M x (0, +cx)). 
Define 

a\ N 'an '1+ P/A4 



mo. 
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Then ([23|-([2l]) imply 

/,,(a(i?)) > 0, when < R < Ma^; 

f^{a{R))<0, when R>Ma„. 

Besides that, one has P^{t,x,R) > 0, thus fjj{a{R))P^{t,x, R) will change sign 
for different R. On the other hand, when 77 <C 1, ( 3.15^ , ( |3.f 6| imply that for 
Vi? e (0, +00), fr,{a{R))P^{t,x, R) has the same sign as /^(a(0))P='=(t, x, 0), which 
is a contradiction. Therefore, /,,(a(0))P* (t, a;, 0) — > 0, and as 77 — 0, 

(3.18) fn{a{R))P^{t,x,R) ^0, Vi?e(0,+oo). 

Then the definition of /(a) in ([23])-([Z4]) gives that if i? 7^ Mq, P^{t,x,R) 0, 
which implies when — ^ 0, 

(3.19) P^{x,t,m) = P^S{m~ Mo). 

4. Keller-Segel limit and connections to the original PBMFT 



In this section, we derive the Keller-Segel limit from (3.5)-(3.8) by taking into 
account the different physical time scales of the tumbling, adaptation and experi- 
mental observations. We shall also connect the new moment system to the original 
PBMFT developed in 25 . Moreover, a numerical comparison of the moment sys- 
tem (3.5|-(3.8| with SPECS is provided in the environment of spatial-temporally 
varying concentration. 

4.1. Keller-Segal limit by the parabolic scaling. We nondimensionalize the 
moment system (3.5)-(3.8) by letting 

t — Tt, X — Lx, vq — sqVq, 

where T, L are temporal and spatial scales of the system respectively. Then 

Jp — SQt/p, Jq — SqJqj 

and the system becomes (after dropping the ) 
1 dp sq dJp 

So dJp _ 2 dp si So sodZ 
T dt ^ '"dxL T,-^" T.dm^^'^ ^^'^"^^ 
Idq 

rm 

sodJq__ 2dqsl So OF 

T dt ~ ''''dxL'^ ''T2 dm'"" ^■^"'"T2 ""Ti dm 
where Ti, T2 are the average run and adaptation time scales respectively. 

For E. coll., the average run time is at the order of Is, the adaptation time 
is approximately 10s ~ 100s, and according to the experiment in 34 , the system 
time scale when Keller-Segel equation is valid is about 1000s. Therefore, we consider 



Sq dJq ^_ 
L dx T2 



Fp, 



So 



So dZ 



M{Jg - MJ, 



So 
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the long time regime, where the tumbhng frequency becomes large (the so-called 
parabohc scahng). Let 



(4.1) 



Ti 

i/so 



L/sq 



1, and 



T _ 1 
L/so £ 



Then (3.5l-(3.8) become 



(4.2) e 



dp 

dt 



Jo 



djp 

dx 

_ ^dp Z 



1 f)V 



dF Z 1 dZ 

— {J, ~ MJp) - -Jg - -^M{Jg - MJp). 
dm E E am 



Consider the following asymptotic expansion 



M = Mo +eMi + 
Z = Zo + eZi + -- 



Jp — ^ ~(~" ^J^p ^ H~ ■ ■ ■ 

Jq — (/g ^ ^Jq ^ • ■ ■ 

F = i^o + e-Fi + • • ■ 



Matching the 0(l/e) terms in (4.3) and (4.5) gives 



^oJr = |^(A^oJr-^r)' and ZoJf = |^Mo(Moj(°)-jf ), 



which implies 



Mo4°) = jf, and j(") = j(o)^o. 



Hence the 0(1) term in (4.4 1 yields 



= 0. 



Equating the 0(e) terms in (4.2) and ( |4.4[ ) produces 



(4.6) 



5pW 
9i 



— and 
ox 



OJq (0) 
dt --^^^'P ■ 



Putting together the 0(1) terms in (4.3) and (4.5) brings 



(4.7) 
(4.8) 



-ZoJ«-|^Afo(j(^)-MoJ«) = 0. 
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The above two equations imply 

dx 



(4.9) 



7-^1? 



dx 
14 



dm 



•Jq 



dx 



dx / 



i9m y (9x 



Z, 



-2„,2 



dm dx 



By (3.3), 



Therefore 



(4.10) 



,(0) _ 



(0) 



'^0 ^0 



9pW 
dx 



'7-2„.2 
^0 



Mo 



dx ' 



Substituting ( |4.10p into (4.6 1 gives the K-S equation 



(4.11) 



djfl 
dt 



'dx 



Zn 



dp^\ _ 
~d^ ^ 



'dx 



, dZp dMp 
dm dx ' 



,(0) 



Using (2.6), ( 3.17[ ) and Mq — Mao, Zq ~ 2;(Ma„), the K-S equation becomes 



(4.12) 
with xo 



dt 



dx 



dx 



[XoP 



(0) 



dfo. 
dx ' 



ao). 



Remark. 1. Instead of (4.1), if we consider 

L/so ' L/so 
then the rescaled system becomes 



and 



T 



d_p 
dt 

' dt 

dq 

'di 

dJ, 



djp 

dx 



dx 



-Jo 



1 dZ 

e dm 



(J, - MJp), 



djq 

dx 



.dq 



1 



-Fp, 



1 dF Z 1 dZ 

-FJ„ + —T^iJq - MJp) - - J, - -T-M{Jq 

ne dm e e dm 



dt dx ne 

When K < 0(1/2), carrying on similar asymptotic expansion will produce the same 
Keller-Segel limit (4.12) as e 0. This indicates that when the adaptation time is 
shorter than ^/TTi, the Keller-Segel equation is valid for E. coli chemotaxis. 

2. The velocity scale of individual bacteria is sq- The temporal and spacial scales 
of the system we consider are T and L respectively, therefore the velocity scale of 
the drift velocity Vd — Jp/p is L/T. The equation (4.1) implies Vd/sQ ^ 0{e), 
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which means that in the regime where K-S equation is vahd, the drift velocity is 
much smaller than the moving velocity of individual bacteria. 

4.2. Connection to the original PBMFT. We shall show that, under certain 
assumptions, the moment system (3.5|-(3.8| gets to the original PBMFT in [25|. 



Especially, one of the equations delivers the important physical assumption eqn. 



(3) in 25 



The macroscopic quantities in the PBMFT in [25] are the cell densities to the 
right P+ and to the left P~ , and the total density ps — + P~ and cell flux 
Js — vo{P^ — P^); the average methylation level to the right Af+ and to the left 
M", the methylation difference AM^ = i(Af+ - Af") and the average methylation 

_ M+P++AI- P- ' ' 



Ms = 

(4.13) 
(4.14) 
(4.15) 



The model in 



25 



is 



dps 

dt 
dJs^ 

dt 

OMs 
dt 



dJs^ 
dx ' 
2dps 



'dx 
J, dM. 



ZJs 



dx 



~ vo-—AMsPs, 
om 

— — (wqAA/^Ps 
Ps dx 



together with the physical assumption 



(4.16) 



AM, 



dM, 
dx 



which physically means AMs is approximated by the methylation level difference 
in the mean methylation field Ms{x,t) over the average run length voZ~^, due to 
the fact that the direction of motion is randomized during each tumble event. 



We firstly discuss the connections of (3.1 1-(3.2 ) to the macroscopic quantities in 



(4.13 )-(4.16 1. By definition, one has 



(4.17) 



P+ + P- 



Ps 



Jp=Vo{P^ 



p- 



Js 



(4.18) q = M+P+ + M-p- = MsPs, 

(4.19) Jq = vo{M+P+ - M-p-) = MsJs 



vg^Ms I 2 

Ps ^ 



Assumption B. Defining the drift velocity Vd — Js/Ps, we assume 

\vd\ < "0- 

This is a fair assumption since in experiments an individual cell usually travels 
at a much higher speed than that of chcmotaxis. 



Applying Assumption [B] in (4.19) gives 
(4.20) 



Jg-MJp 



VqAMsPs 
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by which, (3.5), (3.6) and (3.7) become (4.13), (4.14) and (4.15) respectively. In 



particular, substituting (4.20) into (3.8) produces 



(4.21) ^(M,Js+voAMsp,) 
at 



2 d 



{psMs 



FJs + va^AMsPs) - (msZJs + fo^^^^AM.p, 



Using ( |4.13[ )-( |4.15[ ) yields 
(4.22) vqPs 



dt ps dx ps dx 
+ vqAAL 



(vQpAMs) 



dx 



-oPs^+voPsAMJ—-Z 



Assumption [B] formally implies the terms containing Jg are relatively small in 



(4.22 ), and by a quasi-static approximation dAMs/dt « 0, one has the sum of last 



three terms is approximately zero in (4.22), and thus 

dM, 1 



AAL 



dx + 

Ps ax 



OF 
dm 



- z 



When Z » 



the above equation leads to the important physical 



assumption (4.16) 
(4.23) 



AM 



dM 
dx 



Z ^wo, 



which recovers the PBMFT model in 25 



Remark. The model (4.14)-(4.16) is a nonlinear advection-diffusion system. The 



new moment system (3.5)-(3.8) evolves only linear advection terms and nonlinear 



reactions, and thus the numerical methods for such a system is well studied 22 



4.3. Numerical comparison to SPECS. To show the validity of the moment 



system (3.5)-(3.8), numerical comparisons to SPECS will be presented in this sub- 



section. We choose spatial-temporal varying environment to show how the intra- 
cellular dynamics affects the E. coli behaviors at the population level. Specifically 
it presents a pattern of traveling attractant concentration wave, in which an inter- 
esting reversal of chemotaxis group velocity was revealed in ^ . 

We consider a circular channel with the travelling wave concentration given by 
[L]{x, t) — [L]o + [L\a + sin[^(a; — ut)]. The wavelength A is fixed to be the length 
of the channel, while the wave velocity u can be tuned. The steady state profiles 



of all the macroscopic quantities in (3.5)-(3.8) and corresponding SPECS results 



are compared in Figure [3] The results from SPECS and the moment system are 
quantitatively consistent. It can be noticed that, when the concentration changes 
slowly [u = 0.4/im/s), the profile of M can catch up with the target value Ma^ 
(defined by a(\L\, Mag) = a^), while in the fast- varying environment {u = 8fj,m/s) 
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200 400 




Figure 3. Numerical comparison between the new moment sys- 
tem of PBMFT and SPECS. The steady state profiles of p: (a, 
e), Jp-. (b, f), M — q/p: (c, g), Jq. (d, h) when the traveling 
wave speeds are u — 8fim/s and u — OApm/s respectively. In the 
subfigures, red lines and dots are from SPECS (red lines in b, d, 
f, h are the smoothed results of the red dots), while blue lines are 
from the new moment system of PBMFT. Parameters used here 
are [L]q = 500/iM, [L]a = WOpM, A = 800^*™. 20000 cells are 
simulated in SPECS. 

there is a lag in phase between M and Mag ■ This difference is caused by the slow 
adaptation rate of cell and it also leads to the difference in the profiles of p and even 
chemotaxis velocity; we refer interested readers to [25' for more detailed discussions 
and physical explanations. 



5. Two DIMENSIONAL MEAN-FIELD MODEL 

In this section, we derive the two-dimensional moment system of PBMFT based 
on a formal argument using the point-mass assumption in mcthylation and the 
minimization principle proposed in [18| . 

In two dimensions, v — z;o(cos^, sin0), where vq is the velocity magnitude. 



P{t, X, V, m) in (2.7) can be rewritten as P{t, x, 9, m). z{m, [L], 6, 6') is the tumbling 



rate from 9' to 9. The tumbling term Q{P,z) in (2.8) becomes 
(5.1) 

Q{P,z)= 4- z{m,[L],9,9')P{t,x,9\m)d9' - -f z{m,[L],9' ,9)d9'P{t,x,9,m), 
Jv Jv 

where V — [0, 27r) and j-= ^ Jy. According to (2.6|, z{m, [L], 9, 9') is independent 
of 9 and thus we denote it by z(m, [L]). 
Define 

(5.2) g{t,x,9)^ J P{t,x,9,m)dm, h{t,x,9) = J mP{t,x,9,m) dm; 
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(5.3) 



M{t, X, 9) = 



h{t,x,9) 



M{t, x) = 



Jy h{t, X, 9) d9 



g{t,x,9)' ' \yg{t,x,e)Ae' 

and the density, density flux, momentum (in m), and momentum flux as follows 

(5.4) p{t, x)^ -f g{t, X, 9) d0, Jp(t, x)^ 4- vg{t, x, 9) d9- 

Jv Jv 

(5.5) q{t,x) = -f h(t,x,9)d9, Jq{t,x) ^ -f vh{t,x,9)d9. 

Jv Jv 

We assume 

(5.6) P{t,x,9,m) = g{t,x,9)S{m~ M{t,x,9)). 



This assumption is motivated by (3.19) in one dimension, which could be formally 
understood as the limit of ku — >■ +oo. 
Denote 

3^7 a„ ay 



Integrating (2.7) with respect to m yields 



(5.7) dtg = -v-V^g+ / z{Mi9'),[L])g{t,x,9') d9' - z{M{9),[L])g{t,x,9). 

Jv 

Integrating (5.7) with respect to 9 gives the equation for density, 

dp{t, x) 



(5.8) 



dt 



dJo 



Multiplying (5.7) by v and integrating with respect to 9 produce 

, )d9 

dZ 



- 'j v®vV^gd9~ j vz{M{9),[L\)g{x,t,9) 

(5.9) --J v®v\I^gd9- 

Jv 



v{z{M, [L]) + -T—{M{9) - M))q{x,t,9) d9 
am 



dZ 



= -f v(E)vWa,gd9 - ZJp- ^{Jg- MJp), 
Jv dm ^ ^ 

where we have used the first order Taylor expansion in the second step. 
Integrating mx(2.7) with respect to m brings 

(5.10) dth^ -v-VJi + f{M{9),[L])g{9)+ i z{M{9'),[L\)g{9')M{9')d9' 

Jv 

-z{M{9),[L])g{9)M{9). 



Integrating (5.10) with respect to 6*, and using the definition in (5.3) give 
dq{x,t) 

= — V K • Jq -I- ■ 

iv 



dt 



(5.11) 



4 f{M{9),[L])g{x,t,9)d9 
Jv 

£ {f{M, [L]) + ^{M{9) - M))g(x, t, 9) d9 



-Va, • Jq + Fp. 
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Finally, we integrate t;x(5.10) with respect to 9, 
(5.12) 

dJa 



'^ = -j^v®vVJide+ j^v{f{M{e), [L]) - z{Mid), [L])Mi9))g{d) de 

r 

= - v(S)v- v^hde + Fjp + Q:^^{Jq- MJp) - zjq - - ^^-^p)- 

In order to close the system, especially to find equations for p, J, q and Jg, we 
need to a constitutive relation that represents v®v-V ^gdO and v®v-^ ^liAB 
by p, J, q and J^. The minimization principle proposed in 18 (mathematically a 
projection of g and h on the linear space spanned by 1 and v) gives 

git,x,e) gi{t,x) + gc{t,x)cos9 + gsit,x)sme, 
h{t, X, 0) ss /ii(<, x) + /ic(^j a') cos + /is(t, a;) sin 6*. 



Then from (3.11, (3.2), 



p{t,x) « 




Jp{t,x) « 


/" 


(7(t,a;) w 


/(' 

Jv 


Jq{t,x) « 


/" 



/ii + he cos 9 + hg sin d6' = /ii 



-f f (/ii + hf.cos6 + hg suiO) d9 = —{he, hg)'^ . 
Jv 2 



Therefore, expressing gi, g^, gg, Mi, M^, Mg by p, Jp, q, Jg, we find 



91 = P, 
hi = q. 



9c = 

hr = 





9s = 
















'^Jq.y 




Vo 



where we denote x = (x, y) and Jx and Jy are the x and ?/ components of J. Hence, 



(5.14) 
(5.15) 

Furthermore, since 
(5.16) 



r 2 2 

/ v(8)v- V/id6'«^V/ii = ^Vg. 

2 2 



M 



we are able to close the system (5.8|, (5.9), (5.111, (5.12) using (5.13). 

In summary, ( |5ll| , ( [5l2| together with dsH] ), dsls] ) and 

give us a two-dimensional moment system of PBMFT that is similar to (3.5|-(3.8|. 
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6. Discussion and conclusion 



In seek a model at the population level that incorporates intracellular pathway 
dynamics, we build a new moment system of PBMFT in this paper by using moment 
closure technique in kinetic theory under the assumption that the methylation level 
is locally concentrated. The new system is hyperbolic with linear convection terms. 
Under certain assumptions on the drift velocity, the new system recovers the original 
model in 25 . Especially the assumption on the methylation difference made in 



25 can be understood explicitly in this moment system. We show that when the 
average run time is much shorter than that of the population dynamics (parabolic 
scaling), the hydrodynamic limit of the moment system can be described by the 
Keller-Segal model. We also present numerical evidence to show the quantitative 
agreement of the moment system with SPECS ([l9]). 

We remark that the idea of incorporating the underlying signaling dynamics into 
the classical population level chemotaxi description has appeared in the pioneer 



works of Othmer et al 12 13 33 . Here, the internal dynamics follows the physical 



model proposed in 25 , which results in a closure strategy different from [12[|13 33 



The major differences between the kinetic model in 112 13,33 and the one used 



here are: the methylation rate function is nonlinear (in the methylation level) in 



25 while linear in 12 13 33 ; the tumbling frequency (2.6) follows the results of 



recent physical studies on chemotaxis ( 21 23]). 

Another interesting behavior related to the chemo-sensory system of bacteria 
is the "volcano effect" observed numerically in [oj. It may be also important to 
study this phenomena in a more physical way and understand the communications 
between bacteria using the moment closure technique introduced in this paper, 
which will be our future study. 



Appendix 



We give a systematic way of obtaining systems with higher order moments by 
further introducing 
(6.1) 

e{x,t)= {m- Mf{P+ + p-)dm, Je{x,t) ^ vq {m- Mf{P+ - p-)dm, 



and finding a system of six variables p, g, e, Jp, Jq and Jg. The calculations are 



almost the same as those of deriving (3.51, (3.6), (3.7) and (3 



Define M+, M", M, P±, Z, F and |^ the same as in (|33|, (|0), and 
introduce 



dm'^ ' 



g2p ^ g2f 

dm? dw? '™=A^ 
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The equation for density p is the same as in (3.5 ). Since we keep more terms in the 
Taylor approximation for z(rn), the equation for the density flux Jp becomes 



5Jp 
dt 



'dp 



2dp 
(p+ - p-)dm 

2^9 



j z{m){P+ - p-)dm 



z{M) 



dz 
dm 



{m-M) + -—^ 
m=M 2 am^ 



m=M 



(m - M)' 



dZ dz Id^Z 



Similarly, the equation for q is 

dq 
dt ' 



9a; 2 dm^ 



and for Jg is 



dt 



3W f)7 
FJp + ^{Jg- MJp) - ZJg - ^M{Jg - MJp) 



d^F jyjd^ 2 — ] 
dm^ dm^ dm ) 



J. 



We need two additional equations for e and Jg. Since 

e{x, t) = J - M f{P+ + p-)dm = J rn^{P+ + P^) dm - 2Mq + 
= J m^iP+ + P^)dm- Mq, 



Jeix,t)^vo {m-MY{P+ -P-)Am^va / m^(P+ - P") dm - 2Af J, + MVp, 



we can get the equation for e by multiplying both sides of (2.1)+(2.1) by m^ and 
integrating with respect to m, 



d{e + Mq) d{Je 


+ 2MJg 


- M\jp) 


dt 


dx 




d{Je 


+ 2MJq 


- M^Jp) 




dx 




d{J, 


+ 2MJq 


- M^Jp) 




dx 




Id^ 




(m — 


^ 2 


dm^ 


d{Je 


+ 2MJq 


- M\jp) 



, d{f{a){P+ + P-)) 
dm 

+ 2 / m/(a)(P+ + P")dm 
d{mf{a)) 



dm 



MF 



dm 



=M 



(m - M) 



dx 



{m~MY \{P+ +p-)dm 



8'^ F 8F 
2MFp+(Mj4 + 2|^)e. 

am^ dm 
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The equation for Jg can be obtained by multiplying both sides of (2.1l-(2.1) by 
Vqiti^ and integrating with respect to m, 

d{Je + 2MJq - APjp) 



dt 

^d{e + Mq) 
"° dx 

^d{e + Mq) 
"° dx 

^d{e + Mq) 
"° 



, d{f{a){P+ - P-))) ^ f 

vq m — ^ — am ~ vq z[m)m iF — F ) am 

am J 

+ 2wo / TO/(a)(P+-P")dm-Wo / z{m)m^{P+ - P~) Am 



+ 2vq / (mJ{a))\m=M + 



d{mj) 



=M 



(to - M) 



2 dm^ 
d{z(m)m^) 



dm 



(to - MY (P+ - P )dm-vo / {z{m)m 



V, 



dm 

^d{e + Mq) 



m=M 

(TO-M) + i^'("("^)"^') 

m=M 2 



=M 



dm^ 



m=M 



{m- My ]{P+ - p-)dm 



dF 



0' 



dx 



+ 2MFJp + 2{M— + F) (J, - MJp) + {M 



d^F dF. 



dZ 



l.d^Z 



dm^ dr 
dZ 



M'ZJp {2MZ + M^—) (J, - MJp) ^ ^ij^M' + 4M— + 2Z) J, 



2 dm'^ 



dm 



All these six equations for p, Jp, q, Jq, e, Je together give us a closed moment 
system. 
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